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We extend the formalism of the statistical theory developed for the 2D Euler equation to 
the case of shallow water system. Relaxation equations towards the maximum entropy state 
are proposed, which provide a parametrization of sub-grid scale eddies in 2D compressible 
turbulence. 

1 Introduction 

Two-dimensional flows with high Reynolds numbers have the striking property of organizing 
spontaneously into large scale coherent vortices [Q]. These vortices are common features of 
geophysical and astrophysical flows with the well-known example of Jupiter's Great Red 
Spot, a huge vortex persisting for more than three centuries in a turbulent shear between 
two zonal jets. These vortices share some common features with stellar systems like elliptical 
or spherical galaxies that form after a phase of "violent relaxation" [§, ||, f|]. They can also 
have applications in the process of planet formation which may have begun inside persistent 
gaseous vortices born out of the protoplanetary nebula|5], |(| [5], ||]. Understanding and 
predicting the structure of these organized states is still a challenging problem. 

Since the dynamics of these systems is highly nonlinear, a deterministic description of 
the flow for late times is impossible and one must recourse to statistical methods. The 
statistical mechanics of 2D flows was first considered by |J Onsager (1949), followed by [fTo|l 
Joyce &; Montgomery (1973), in the point vortex approximation. The more realistic case 



of continuous vorticity fields was later on treated by |ll[ Kuzmin (1982), [|12[ Miller (1990) 
and (l3|] Robert & Sommeria (1991). 

The statistical theory of the Euler equation provides a systematic framework to tackle 
the problem of self-organization in 2D flows. At a microscopic level, complex, inviscid, 
deformation of the vorticity field creates an intricate filamentation; however, if we introduce 
a macroscopic level of description (a "coarse-graining") it can be shown that an overwhelming 
majority of these microscopic configurations are close to a macroscopic state (the statistical 
equilibrium or Gibbs state) obtained by maximizing a "mixing entropy" while accounting 
for all the constraints of the Euler equation (the conservation of energy and of the detailed 
distribution of vorticity). The maximum entropy theory generally predicts a nonlinear 
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relationship between vorticity and stream function which respects the properties of the 
inviscid dynamics. 

The predictions of the statistical theory [y = 0) have been tested in a large number of 
numerical simulations and fluid laboratory experiments at high Reynolds numbers [y — > 0). 
Good agreement is obtained in many cases, as long as stirring is sufficiently intense to lead 
to equilibrium before significant influence of viscous effects H, pTof] . 

However the Euler equation itself is of limited scope for applications to atmospheric or 
oceanic motion, where the Coriolis force and density stratification have a strong influence. A 
first step in this direction is provided by the quasi-geostrophic model. The application of the 
statistical theory to this case is straightforward as the flow is still assumed non-divergent, 
and the vorticity is just replaced by a potential vorticity. This has been discussed in several 



papers @, [TJ, |J, [lj . 

We here extend the theory to the more general case of shallow water system, a "compress- 
ible" 2D flow, whose properties are recalled in section ||[ Several numerical computations, 
for instance |^, indicate inertial organisation of vortical motion into coherent vor- 

tices, like with the incompressible Euler equations. In these flows dominated by vortical 
motion, the influence of density waves is weak, and the free surface tends to be controled by 
the vortical motion through a balance condition generalizing the quasi-geostrophic balance. 

In section [3] the procedure of [13] Robert and Sommeria (1991) is extended to the 
shallow water system with discussion of simplified cases in section ||. We still assume 
that the vorticity field creates intricate filamentation but the divergence and water height 
(surface density) fields are still smooth. These assumptions are justified for flows dominated 
by vortical motion at moderate Mach numbers, for which the generation of shocks is not 
effective. The relaxation toward equilibrium is discussed in section |E| providing practical 
methods for determining the statistical equilibrium as well as sub-grid scale modeling of 
turbulence in a shallow water system. Finally the case of particular geometries is discussed 
in section @. 



2 The shallow water equations 

Consider a fluid layer with thickness h(x, y, t) submitted to a gravity g on a rotating planet. 
We assume that the layer is thin with respect to the characteristic length scale of the 
horizontal motion. In that case, the velocity field u(x, y, t) can be assumed two-dimensional 
and the vorticity w = we z = V A u is directed along the vertical axis. We shall assume for 
simplicity a plane geometry, with rotation vector Q directed along the vertical, but extension 
to a spherical geometry would be straigthforward (we introduce the Coriolis effect but no 
centrifugal force as the latter is incorporate in the gravity of the planet). The time evolution 
of these quantities is determined by the shallow water equations (see, e.g., j23| ): 

|£ + V-(fcu) = (1) 

— + (w+2Q) Au = -VB (2) 

Here the usual advective term u.Vu has been expressed using the well-known identity of 
vector analysis u.Vu = V(u 2 /2) +wA u, and the term u 2 /2 incorporated in the Bernouilli 
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function 



u 2 



B = gh + —, (3) 

together with the pressure term gh. Note that the shallow water system can be viewed as a 
2D flow of a compressible gas with density h and equation of state p = gh 2 /2, and our results 
could be readily generalized to 2D compressible adiabatic flows. We shall often refer to the 
shallow water system as the "compressible case", by opposition with the "incompressible" 
case, for which (|j) is replaced by V.u = . 

One can easily check that the potential vorticity (PV) 

q = ^r~ (4) 

is conserved for each fluid parcel, i.e: 

$-!-V.-0 ,5, 

Each mass element hd 2 r is conserved during the course of the evolution. Together with (||), 
this implies the conservation of the Casimirs 



C f = J f(q)hd z v (6) 

where / is any continuous function of the potential vorticity. In particular, the moments 
r n of the potential vorticity are conserved 



T n = J q n hd 2 r (7) 



The moments n = 0,1,2 are respectively the total mass M, the circulation V and the PV 
enstrophy T2. The energy 

E = J h ^d 2 v + \J gh 2 d 2 v (8) 

involving a kinetic and potential part, is also a conserved quantity. Note finally that for a 
multiply connected domain, like the annulus or channel discussed in section ||, the circulation 
along each boundary is conserved, in addition to V . 

It will be convenient in the sequel to use a Helmholtz decomposition of the momentum 
hu into a purely rotational and a purely divergent part 

hu = -e z A Vip + Vcp (9) 

where tp and are defined as solutions of the Poisson equations 

= —V A (hu), tp = const, on each boundary (10) 

A(p = —V • (hu), d(p/dC, = on each boundary (11) 

where the conditions at the domain boundary (with normal coordinate C ) are the conse- 
quences of the impermeability condition. We here consider a domain with a single (outer) 
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boundary, so we can take i[> = at this boundary (as if) is denned within an arbitrary gauge 
constant). 

For a steady solution, the mass conservation (|]) reduces to V • (hu) = 0, so that 4> = 

and 

hu = —e z A Vtp (steady), (12) 

and equation (||) reduces to u • Vq = 0, which implies that q is a function F of the stream 
function ip 

q = Fty) (13) 

Finally, equation (§) reduces to 

(u+20)Au = -V_B (14) 

Taking the dot product with u, we obtain u • VB = or, equivalently, B = B(ip). This is 
known as the Bernouilli theorem. Substituting for (12) in equation (|il|), we obtain a specific 



relationship between the potential vorticity q and the Bernouilli function B in the form: 

Small perturbations to a state of rest, with uniform thickness H , satisfy linearized 
equations with two branches of solutions. For small scales, these are the usual surface 
waves on one hand, with purely divergent velocity and propagation speed c = (gH) 1 ^ 2 , and 
steady vortical divergenceless modes on the other hand. In nonlinear regimes, these two 
modes interact. Vortical motion with scale / and typical vorticity u fluctuates on time scale 
tu^ 1 , so it emits waves at wavelength A ~ c/u>, hence X/l is the inverse of the Mach number 
c/u based on the local velocity u ~ ujI. Therefore we expect that for the case of small Mach 
numbers that we shall consider, velocity divergence and free surface deformation are much 
smoother than the vorticity field (their wavelength is much larger). 



3 The maximum entropy theory 
3.1 General principles and notations 

For slow velocities, the shallow water system reduces to the quasi-geostrophic (QG) equa- 
tions, with h ~ cte , such that (|]) reduces to the incompressibility condition V.u = . Then 
the velocity field remains regular for any time, but it generally develops complex fine scale 
vorticity filaments so that a deterministic description of the flow would require a rapidly 
increasing amount of information as time goes on. The idea of the statistical theory is to 
give up such a deterministic description and refer to a probabilistic description. Therefore, 
the exact knowledge of the "fine-grained", or microscopic potential vorticity field is replaced 
by the probability density (area fraction) p(r, a) of finding the potential vorticity level a at 
position r. 

For the more general shallow water system, the inviscid dynamics generally leads to 
singularities (shocks), with associated energy dissipation (even in the absence of viscos- 
ity). This is a source of fundamental mathematical difficulty for the generalization of the 
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equilibrium statistical mechanics initially developed for the Euler equations or QG system. 
Nevertheless, for the case of small Mach numbers that we consider, shocks occur only after 
a very long time (due to non-linear steepening of surface waves), and they involve a weak 
energy dissipation, since most of the energy remains in the vortical motion. Furthermore 
fine scale vorticity fluctuations behave like in the QG system, and only interact with surface 
waves and flow divergence at much larger scale, as discussed above. We shall therefore as- 
sume that vorticity fluctuates at small scale, but divergence is smooth as well as the height 
h. 

We still describe the local vorticity fluctuations by the probability p(r, a) of finding the 
potential vorticity level a in a small neighborhood of the position r. The normalization 
condition yields at each point 

' p(r,a)da = l (16) 



The locally averaged field of potential vorticity is expressed in terms of the probability 
density in the form 

q = J p(r,a)ada (17) 

This locally averaged field is called the macroscopic, or coarse-grained, potential vorticity. 
A macroscopic state is fully defined by p(r, cr), the height field h(r) and the flow divergence 
(assumed without microscopic fluctuations). The velocity field u(r) can be also considered 
as smooth, as it integrates the vorticity fluctuations, and can be deduced by integration 
from its divergence and vorticity U = qh — 2Q . The energy (||) depends only on this smooth 
field, with negligible influence of local fluctuations, like in the incompressible case 

The conservation of the Casimirs @ is equivalent to the conservation of the global 
distribution of potential vorticity (i.e the total area of each level of potential vorticity 
ponderated by h): 

7 ((t) = / p(r,a)h(r)d 2 r (18) 



The microscopic moments of potential vorticity can be written: 

r n = J ~f{a)a n da = J q^h(r)d 2 r (19) 



where 



q n = J p(r,a)a n da (20) 

and are conserved during an inviscid evolution. Note that for n > 2, the macroscopic 
moments of potential vorticity F^ 9 ' = J q n h(r)d 2 r are not conserved, as there are partly 
transfered into fine-grained fluctuations. 
The mixing entropy 

S = - J p(r, a) lnp(r, a)h(r)d 2 rda (21) 

measures the number of microscopic configurations associated with the same macroscopic 
field of potential vorticity. The dependence in p is the same as for the incompressible 



case [13], and the factor h(r) is introduced to insure that entropy is conserved by mere 
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macroscopic displacement of fluid parcels. Indeed the mass element h(r)d 2 v is conserved 
by fluid particle displacement instead of the surface element d 2 r in the incompressible case. 
At equilibrium, the system is expected to be in the most probable (i.e. most mixed) state 
consistent with the constraints of the Euler equation. The entropy ( |2l| ) has been justified 
by rigorous mathematical arguments (in the incompressible case) but other forms have been 
recently proposed (see discussion in |l^]). Therefore, we shall consider a general form of 
entropy 

S = - J s(p{r,a))h(r)d 2 rda (22) 
and find that equation (|2l|) is the only expression leading to a well defined entropy extremum. 

3.2 First order variations 

According to the previous discussion, the most probable macroscopic state is obtained by 



maximizing the entropy (|22[) with fixed energy (g), global vorticity distribution ((18|) and 
local normalization (|l6|). This problem is treated by introducing Lagrange multipliers, so 
that the first variations satisfy 

5S - p5E - J a(a)5-f(a)da - J C(r)5^J p(r, a)do\ hd 2 r = 0. (23) 

The Lagrange multipliers are respectively the "inverse temperature" (3, the "chemical poten- 
tial" a(a) of species a, and C( r )- 

We shall take h, ph and V • u as independent variables characterizing the macroscopic 
state. Then, it is easy to establish, by differentiating respectively (|2l|), ( |l8| ) and (16), that: 

5S = J [ps'(p) - s{p)]5hd 2 rda - J s'{p)5{ph)d 2 rda (24) 

6i(a) = J 5{ph)d 2 r (25) 

hsfj p(r,a)da^j = J 5{ph)da - J p5hda (26) 

The variations of the energy are conveniently expressed in terms of the Bernouilli function 

B, 

5E= [ B5hd 2 T + [ hu ■ 5ud 2 r (27) 



Then, using the Helmholtz decomposition ([|) for the momentum hu , the second integral 
can be rewritten 

J hu ■ 5ud 2 r = - J (W A 5u) ■ zd 2 r + Jvcp- 5ud 2 r (28) 

Integrating by parts with the identities V A (ipdu) = ipV A (<5u) + Vi/j A 5u and V • (<M U ) = 
0V • (5u) + V(p ■ 5u, and using the boundary conditions for ip and (j) , we obtain 

J hu ■ 5ud 2 r = J il)5u5d 2 Y - J (jN ■ (5u)d 2 r (29) 
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Using (H) (|) and (0), we have finally 

5E = J Bp6hd 2 rda + J ^a5(ph)d 2 rda - J 4>5{V ■ u)d 2 r (30) 

The variation (23) vanishes for any small changes of the variables only if the coefficient 
of each independent variable vanishes: 

8{ph) : s'{p) = -p<nj) - a{a) - C(r) (31) 

Sh : s'(p) - S -^- = f3B- C(r) (32) 

<5(V • u) : = (33) 



Note that the right hand side of equation (|32j) is independent of a. This implies that the 
term on the left hand side must be a constant (that we can take equal to 1 without loss of 
generality) :s'(/>) — s(p)/p = 1. This equation is easily integrated in s{p) = Ap + p In p where 
A is an integration constant. When substituted in equation using (|l9|) , this yields 



S = - J p(r, a) In p(r, a)h(r)d 2 rda - AM (34) 



which is just the entropy (|2l|) up to an additive constant term AM (which we can take equal 
to zero without loss of generality). Therefore, the entropy ( |2l| ) is the only functional of the 
form (^2|) for which the maximization problem has a solution. This result is astounding 
because it is obtained without any explicit reference to thermodynamical arguments. 

3.3 The Gibbs states 

Substituting explicitely s(p) = pin/? in equation (31), the optimal probability density can 



be expressed as 

p(r,a) = ^- ) g(*)e-W (35) 

where Z(ip) = e^ r ) +1 and g(a) = e~ a ( a \ The normalization condition (|l~6|) leads to a value 
of the partition function Z of the form 

Z = f g(a)e~^da (36) 



and the locally averaged potential vorticity (|17|) is expressed as a function of ift according 
to: 

J g(a)ae~^da 



J g(a)e~^da 

This can be rewritten 



F(V0 (37) 



1 dlnZ , , 

■TW (38) 



This is the same functional relation as in the case of 2D incompressible Euler flows |13| . 
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Differentiating equation (37) with respect to ip, we check that the variance of the po- 
tential vorticity can be written 



Q2 = T 



or, alternatively (see equation (^H)): 



<12 



1 d 2 lnZ 

P 2 dll> 2 



(39) 



(40) 



Therefore the slope of the function q = F(ip) is directly related to the variance of the 
vorticity distribution. Relation (|39|) has the same form and origin as the "fluctuation- 
dissipation" theorem in statistical field theory, where dq/dip is interpreted as a susceptibility. 
Since qi > 0, we find that the function q = F(ip) is monotonic; it is decreasing for f3 > 
and increasing for f3 < (it is constant for (3 = 0). Another proof of this result is given in 

El- 

Substituting explicitely s'(p) — s(p) / p = 1 in equation (|3^), we have 



D 



InZ 



(41) 



This relation shows that the Bernouilli function plays the role of a free energy in the sta- 
tistical theory. We note that both B and q are functions of ip, while = 0, from (|33|), 
as it should for steady flows. Furthermore, taking the derivative of (|4l|) with respect to ip 
and using (|38|), we check that the relation q = —dB/dtp required for a steady solution of 
the shallow water equation is satisfied. Therefore the flow selected by the purely statisti- 
cal entropy maximization procedure does not evolve anymore by the flow evolution, so the 
statistical theory is indeed consistent with the dynamics. 



4 Properties of the Gibbs states in some particular cases: 
4.1 Particular q — if) relationships: 



The Gibbs states are characterized by the relation (|37|) between q and ift, which is always 
monotonic, as shown in the previous section. It is determined by the conservation laws, 
but only indirectly through the set of Lagrange multipliers (3 and a(a). In practice we 
need to discretize the PV levels, and keeping only two levels, q = ctq and q = o\ is often 
representative of more general cases. Then the probability distribution p just depends on 
a single probability p\ of finding the level o\ (for instance), with the probability 1 — p\ 
of finding the complementary level do , i.e. g{a) is the sum of two Dirac function terms, 
g(a) = g\[X5(a — do) + 5(a — cr%)]. This probability p\ is directly related to the PV average 
by q = V\ G \ + (1 — Pi)°U) or reversely Pi = (q — c )/(cji — do). The expression ( |37| ) for q 
reduces to 
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This relation corresponds to the Fermi-Dirac statistics. The two unknown parameters A and 
(3 are indirectly determined by the conserved quantities. The associated Bernouilli function 
(41) becomes 

B = i In 51 - ctoV + j ln {^ + e^( ff °- CTl ^| (43) 
The problem is also greatly simplified in the alternative case for which g{a) is a Gaussian: 



g(a) = g e ^ . (44) 



Then the local probability distribution (|35|) is also a Gaussian, and the corresponding 
Bernouilli function ( [Sl|) is 

B = i ln[g (2ira 2 ) 1/2 } + \a*&tf - a*ip, (45) 

corresponding to a linear relationship: 

q = -pa 2 ij> + a* (46) 



According to (39) the variance of the potential vorticity is then uniform, with value qi 



cr 2 (more generally, all the even momenta of the Gaussian are related to 02 by (q — q) 2n = 
(2n — l)!!o"2 and the odd momenta cancel). 

This Gaussian local probability distribution is obtained by maximizing the entropy (|2l|), 
reducing the constraints of the global distribution 7(0") to its first moments Tq = M, Ti 
and T2. This will be the true Gibbs state for a particular initial distribution 7(0") with 
higher order moments equal to the global moments of this simplified Gibbs state. A linear 
relationship between q and tp can also be obtained for any distribution 7(0") in the limit of 



strong mixing where /3cn/" *C 1, so that (37) can be linearized, as discussed in |24J] Chavanis 
& Sommeria (1996). 

4.2 Unidirectional solutions 

We consider here unidirectional solutions such that u = u(y)e x . The equilibrium relation 
q = —dB/dtp then yields, multiplying each member by hu = dijj/dy, 

^dh + 2Q dip _ Q 
dy h dy 

This condition of geostrophic balance can be readily integrated in 



/ 4Q 

h = H Y~W i) (48) 

A second relation is provided by the expression of the Bernouilli function yielding 

1 / dib \ 2 

1 > =B{i>)-gh (49) 



2ti 2 V dy 
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Combined with the expression ( fig] ) of h and the Gibbs state expression ([H|) of B(ip), 
this yields a first order ordinary differential equation for ip. This equation depends on the 
constant H and Lagrange parameters, for instance g±, j3 and A in the case with two PV 
levels. The constraints on total mass M, global mass of PV level o\ and total energy E, 
indirectly determine these parameters. Note that in reality this solution must be viewed 
in a x-wise translating frame of reference, as discussed in section pi due to the additional 
conservation law for momentum. 



4.3 Axisymmetric solutions 

For axisymmetric solutions u = ug{r)eg where (r, 9) are polar coordinates, and hug = 
—dtp I dr. Then Equation ( |47| ) is replaced by the cyclostrophic balance 

gh dh _ 1 f d ^\ 2 2 n— (50) 

dr hr\dr J dr 

When hug = —dip/dr > (cyclone), h is an increasing function of r, so the vortex core 
is a trough. In the opposite case ug < (anticyclone) , the vortex core is a bump in 
geostrophic regimes. However for large velocities (Rossby number larger than unity), the 

u 2 

term dominates, leading always to a trough. 

The expression of the Bernouilli constant gives again the form (|49|), just replacing y by 
the radius r . Combining this relation with (|50|), one gets a couple of first order ordinary 
differential equations in the variables ifi and h. As in the unidirectional case, the solution 
depends on two constants of integration and the Lagrange parameters, which are g%, f3 and 
A in the case with two PV levels. This solution must be viewed in general in a rotating 
frame of reference, due to the additional conservation of angular momentum , as discussed 
in section 0. 



5 Relaxation equations 

5.1 The Maximum Entropy Production Principle 

Relaxation methods are convenient to compute the statistical equilibrium resulting from 
any initial condition. The aim is to increase the entropy in successive steps while keeping 



constant all the conserved quantities. |25] Turkington & Whitaker (1996) have implemented 
a relaxation method to calculate the Gibbs states obtained with the Euler equations. 
Robert and Sommeria (1992) had previously proposed relaxation equations in the form of 
a parameterization of sub-grid scale eddies which drives the system toward statistical equi- 
librium by a continuous time evolution. Such relaxation equations can be used both as a 
realistic coarse resolution model of the turbulent evolution, and as a method of determina- 
tion of the statistical equilibrium resulting from this evolution (see Q for a short review). 
We here generalize this approach to the shallow water system. 

We first decompose the vorticity u and velocity u into a mean and fluctuating part, 
namely u> = U + u, u = u + u, keeping h smooth. Taking the local average of the shallow 
water equations (|])(§), we get 

^ + V.(/m) = (51) 
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<9u 

— + (tD + 21) A u = -VB - e z A J w (52) 



u 2 



B = gh + — (53) 

where the current J w = cju represents the correlations of the fine-grained fluctuations. 
Although we have neglected the fluctuation energy u 2 in front of u 2 (as well as fluctuations 
of h), we keep the correlations J w = wu, which represent the PV transport by sub-grid- 
scale eddies. This assumption is justified since, denoting e the typical scale of vorticity 
fluctuations, we have u 2 ~ e 2 uJ 2 and cDu ~ ecZJ 2 3> u 2 (while Q ~ uS). 

We deduce an equation for the evolution of the potential vorticity (H), taking the curl 



of ( p2| ) and using (|5T|), 

^- t (hq) + V.(hqu) = -V.J u! (54) 

This equation can be viewed as a local conservation law for the circulation T = J qhd 2 r. 
We shall determine the unknown current by the thermodynamic principle of Maximum 
Entropy Production (MEP) |26|]. For that purpose, we need to consider not only the locally 
average PV field q, but the whole probability distribution p(r, a, t) now evolving with time 
t. The conservation of the global vorticity distribution 7(0-) = / phd 2 r can be written in 
the local form 

^(M+V.(Vu) = -V.J (55) 

where J(r, a, t) is the (unknown) current associated with the level a of potential vorticity. 
Integrating equation (|55|) over all the PV levels a, using ((0|), and comparing with (^|), we 
find the constraint 

J J(r,a,t)da = (56) 
Multiplying (|55|) by a, integrating over all the PV levels, using (|T^) and comparing with 



(|54j), we get 

J J(r, a, t)ada = J w (57) 

We can express the time variation of the energy E = dE/dt in terms of J, using (|8|) and 
([52]), leading to the energy conservation constraint 

E = [ 3ahu ± d 2 rda = , (58) 



where uj_ = Au . Using (21) and (|5q), we similarly express the rate of entropy production 



as 



S = -J 3V(lap)d^rda . (59) 
According to the Maximum Entropy Production Principle, we determine the current J 



which maximizes the rate of entropy production S respecting the constraints E = 0, ( |5q ) 
and / j^da < C(r,a). The last constraint expresses a bound (unknown) on the value of the 
diffusion current. Convexity arguments justify that this bound is always reached so that the 
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inequality can be replaced by an equality. The corresponding condition on first variations 
can be written at each time t: 



5S- (5{t)5E- J C(r, t)5^J Jda^jcfr- J D-\r,t)5^J ^jdad 2 r = 

and leads to a current of the form 

J = -D(r, t)(Vp + P(t)aphu ± - C(r, t)p) 
The Lagrange multiplier C,{v,t) is determined by the constraint (feH), which leads to 



-D(r,t) 



Vp + P(t)p(a -q)hu± 



(60) 
(61) 
(62) 



This optimal current is similar to the one obtained in ordinary incompressible 2D turbulence 
except that the term hu± now replaces Vip. The impermeability condition imposes that the 
normal component of the velocity and of the current vanishes at the wall. We therefore 
have the boundary conditions 

nTu = (on each boundary) (63) 

n.V/O = — (5{t)p{a — q)hriu± (on each boundary) (64) 

where n is a unit vector normal to the wall. 

The diffusion coefficient D is not determined by the MEPP as it depends on the unknown 
bound C on the current. It can be determined by more systematic procedures inspired from 



kinetic theories of plasma physics like in |£7|, [2^] for the incompressible case. For the purpose 
of reaching the Gibbs state, it can simply be chosen arbitrarily. However, we shall show 
below that D must be positive so as to insure entropy increase. 

The conservation of energy ( |5§| ) at any time determines the evolution of the Lagrange 
multiplier (5{t) according to 

/ DVqhu ± d 2 r 



0(t) 



(65) 



J D{q 2 -q 2 ){hu ± ) 2 d 2 r 

We can now provide an explicit form for the vorticity current to introduce back in 
the shallow water equation (|^). Indeed, using ( |62|) and (|l7|), we find 



-D 



Vq + P(t)(Q 2 -t)hu ± 



Substituting ( p6j ) in equation (|52j), we obtain 



<9u _ 

— + (o + 2J1) A u 

ot 



-VB + D 



AVq-(3(t)(Q 2 



hu 



(66) 



(67) 



Since (3(t) < in relevant situations, the last term in ( |67| ) represents a forcing proportional 
to u which compensates the diffusion caused by the term e z A Vq ~ Au. This additional 
term depends on the local PV variance q 2 — q 2 , related to the probability distribution 
p, and we need to keep track of it by solving the probability transport equations (55) in 
addition to the modified shallow water system (51) and (|67|). This set of equations increases 
the entropy (at an optimal rate), while preserving all the conservation laws of the initial 
inviscid shallow water system. We now check that the steady solutions reached by these 
relaxation equations is indeed the Gibbs state. 
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5.2 Relaxation towards the statistical equilibrium 

The entropy production fl59| ) can be written 

J 

-( 

P 



S = - J -(Vp + Pp(a - q)hu ± )d 2 rda + [5 J J(a - q)hu L d 2 rda (6S 



Using (56) and (58), the second integral is seen to cancel out. Substituting for (p2) in the 



first integral, we have 

S = J —d 2 vda (69) 

which is positive provided that D > , and this is clearly a necessary condition to assure 
entropy increase in all cases. A stationary solution S = is such that J = yielding, 
together with (||), 

V(ln p) + p{p - q)Vip = (70) 
For any reference PV level <7o, it writes 

V(ln Po ) + P(a -q)Vi; = (71) 



Substracting (|7fj) and (|7l|), we obtain Vln(p/po) + P(& — ctq)ViP = 0, which is immediately 
integrated into 

pM = W) 9{a)e ~^ (72) 

where Z _1 (r) = po( r )e^ cr °^ ( ' r ^ and g(a) = e A( - a \ A(a) being a constant of integration. 
Therefore, entropy increases until the Gibbs state ( |35| ) is reached, with f3 = limi_ >00 0(t). 

5.3 Simplified cases: 

In the case of two PV levels do and a\, the transport equation (|55|) for the probability 
pi is equivalent to the transport equation for q (since q = do + pi(ci — o"o)) ; which is 
already obtained from the curl of the shallow water equation (|67|) . Therefore the relaxation 
equations reduce to the modified shallow water system 

Bh 

— + V-(hu) = (73) 

+ qhe, A u = -V(y + gh) + D\e, A V? - /3(t)q 2 hu] (74) 

_ = (VAu)„e, + 20 . a-ft-^-a (75) 

J Dq 2 {e z A u) 2 h z d z r 
n.Vq = —0(t)q2hnu± (on each boundary) (77) 
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n.u = (on each boundary) (78) 

where we have omitted the over-bar on u, and the expression ( f75| ) of q2 = q 2 — Q 2 is 
easily obtained for a probability distribution with two values gq and u\. The numerical 
implementation of this system will lead to the two PV levels Gibbs state. 



Stating c/2 = cte instead of the expression ( |75|) yields the Gaussian Gibbs state with linear 



relationship between q and ip. Then the coefficient q20 used in (74) is directly obtained from 
(|76|) . This is sufficient for the purpose of finding the statistical equilibrium, but more refined 
relaxation models can be used as discussed by [^] Kazantsev et al. (1998) in the context 
of QG models. 

5.4 The incompressible limit: 

The case of ordinary 2D incompressible turbulence is recovered in the limit h — ► 1, q — ► u 
and u = — e z A Vip. The relaxation equation ( |67| ) then becomes 



O -i 

^ + (u.V)u = --Vp + J D(Au-/3(t)w 2 u) (79) 

where we have used the well-known identity of vector analysis Au = V(Vu) — V A (V A u) 
which reduces for a two-dimensional incompressible flow to Au = e z A VcZJ. Equation 
(|79|) is valid even if D is space dependant unlike with a usual viscosity term. In previous 
publications this equation was given only in its vorticity form 



^ + V(*n) = V 



D[ VuJ + (3(t)cu 2 Vip 



(80) 



and the equivalence with (|79j) is not obvious at first sights when D is space dependent. At 
equilibrium, we have from ( fZ9| ) the identity 

Au = (3uj 2 u (81) 

which can be deduced directly from the Gibbs state. Indeed for a stationary solution 



ZJ = F(ifj), the previous identity Au = e 2 A VlJ becomes Au = —F'(ip)u equivalent to (81) 
for a Gibbs state thanks to (p9|) . 

We now account for a deformation of the fluid layer but assume that the elevation with 
respect to the average thickness H is weak, so that 

h = H(l + if) with 7] <C 1 (82) 

To first order the flow is incompressible and equation ([[]) reduces to V.u = 0, or equivalently 
u = — e z A V?/> (there is a factor H with the previous definition (^)). In the quasi-geostrophic 
limit of small Rossby numbers lo <C ^, the momentum equation (||) reduces at zero order 
to the geostrophic balance 

gH „ , gH 2 , \ 

u = -^e z AVij or tp = — — ry (83) 
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The expression for the potential vorticity then reduces to 



( = Hq- 2n~uJ + jz 



R 



with the Rossby radius of deformation 



Lr 



20, 



(84) 



(85) 



The term j-^P in fl84| ) creates a shielding of the interaction between vortices (similar to 
the Debye shielding in plasma physics) on a length scale ~ Lr. In the limit 1/Lr — > 0, 
we recover the 2D incompressible equations. For finite Lr we can extend the statistical 
theory of the 2D Euler equations to the QG case by simply replacing the vorticity u; by the 
potential vorticity C M, FTl ||, 111, 



6 The case of circular domains or channel: 

6.1 Statistical equilibrium 

In a disk, the angular momentum 

L = J h(rAu) z d 2 r (86) 
is conserved. This additional constraint can be accounted for by adding a term (3X5L in 



the first order variation (23). We can write 5L = J 5h(r A u) z <i 2 r + / h(e z A r).<5ud 2 r, and 
the second term can be expressed in terms of 5lJ and 5(V.u) by a Helmholtz decomposition 
of h[e z A r) analogous to (|9|), followed by an integration by part. This is analogous to the 



formulae (Eg) (29) used for expressing the energy variation. We can combine the energy and 



momentum variations by defining 

h[u - \(e z A r)] = -e z A Vip' + V<p' (87) 
instead of (||). Adding the new terms in 

©flU y ields the Gibbs state dH)© for the 

velocity u' = u — A(e z A r) seen in the reference frame rotating at angular velocity A. Ac- 
cordingly, the expression of the Bernouilli function must be modified by a term of centrifugal 

force: we must use B'(ip') = gh + ^ A 2 r 2 instead of B{ijj). We find therefore that the 

Gibbs state (its locally averaged velocity field) is a solution of the shallow water equation 
which is steady in a reference frame rotating at a modified angular velocity ft + A. This 
velocity is indirectly determined by the constraint on angular momentum. Note that the 
result can be readily extended to the shallow water system on the sphere. 

In the case of an annulus, the circulation C_ = — / ugdl around the inner wall is an 
additional conserved quantity (the circulation C+ around the outer wall is also conserved, 
but it is related to other conserved quantities, as C + = V — C_, and the conservation of V is 
already included in the P V conservation) . Furthermore we need in general to set ip = ?p- ^ 
at the inner wall (while we can still set tp = at the outer wall). As a consequence 
a boundary term -ip-5C- now appears in the expression (53) for the energy variation. 
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However we can directly set 8C- = 0, canceling this boundary term, without influence on 
the independent variables hp (determining the locally averaged vorticity u = J ahpda), V.u 
and h. Therefore the only modification with respect to the disk is an additional unknown 
■0- in the definition ([[(]) of ip, determined by the corresponding additional constraint C_ . 

The case of a straight channel can be viewed as the limit of an annulus with small width, 
but it can be convenient to treat it in itself. Let us consider a straight channel between two 
walls at coordinates y = ±L y /2 with periodic boundary conditions along the x-direction 
(with domain length L x ). In the absence of Coriolis force (O = 0), the x-wise momentum 

P = I hu T d 2 r (88) 



is conserved (instead of the angular momentum), as well as the circulation C + = — J u x dx 
along the upper wall (y = L y /2 ). The boundary condition (|l0|) defining ip is replaced 
by ip = P/L x at the upper wall y = L y /2 and ip = (for instance) at the lower wall 
y = —Ly/2. Unlike with angular momentum, we cannot express the variation <5P in terms 
of the variations in the independent variables UJ , V.u , h. However we have now an additional 
freedom in the variational problem, as we can add a uniform x-wise velocity —Ue x (use a 
reference frame with velocity Ue x ) without influence on the independent variables ZJ , V.u , 
h. For any choice of U, we can solve the variational problem with the velocity u' = u — Ue 
and corresponding energy E' = E + MU 2 /2 — PU, upper wall circulation C' + = C + — UL 



X 



This yields a Gibbs state (p5|)(|38|)(|4lD, representing a steady solution of the shallow water 
equation in a reference frame moving with velocity Ue x . Among these states, the ones with 
the right value of the momentum P = J hu x d 2 r will be the actual solutions. Families of 
Gibbs states with the same structure translated in the x-direction are obtained, as discussed 
by |29|] Sommeria et al. (1991) in the incompresible case. 

Finally in the case of an infinite domain, the two components of momentum, as well as 
the angular momentum are conserved. This yields to purely translating or purely rotating 
Gibbs states, as discussed by [^] Chavanis and Sommeria (1997) in the incompressible case. 

6.2 Relaxation equations 



Taking the derivative of (88)(|86|) with respect to time and using equations (|5l|) (|52"|) , we 
obtain the constraints 

P = J hJ wy d 2 r = (89) 

L = - J /iJ^.rcPr = (90) 

These constraints can be included in the variational principle (|6C| ) by introducing ap- 
propriate Lagrange multipliers denoted (3(t)U(t) and (3(t)\(t). Then, the results of sec- 
tion H| are generalized simply be replacing the velocity u by the relative velocity u' = 
u — U(t)e x — \(t)e z A r where the time evolution of U(t) and X(t) is obtained by substi- 
tuting the optimal current (|66|), with the above transformation, in the constraints ( |89| ) and 
(|90|). In the case of a channel we have the additional conserved quantity C+ = — J u x dx 
along the upper wall. Using (p^), we readily find that C + = f J^ydx = as the current J w is 
parallel to the wall, so the circulation along each wall is indeed conserved by the relaxation 
equations. 
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